function [Price, V, iteration] = butterflyUnVol( r, sigmaMin, sigmaMax, S0, K1, K2, T, steps, method, cases, tol, NodeM )
%EUIMPLICIT Summary of this function goes here
%   Butterfly option pricing using uncertain volatility method
%   Constant stepping is used
% method: 0=fully implicit (Default), 1 = Crank-Nicolson
% NodeM: Node multiplier, default is 1
% cases: 'w' = worst case long, 'b' = best case long
% tol: tolerance level, default to 10^-6


% Argument check
if nargin < 9
    method = 0; cases = 'w'; tol = 10^-6; NodeM = 1; end
if nargin < 10
    cases = 'w'; tol = 10^-6; NodeM = 1;end
if nargin < 11
     tol = 10^-6; NodeM = 1; end
if nargin < 12
     NodeM = 1; end

if method == 0
    theta = 0;
else
    theta = 1/2;
end

% Generate S grid
S = S0*[0:0.1:0.4, 0.45:0.05:0.8, ...
    0.82:0.02:0.9, 0.91:0.01:1.1, ...
    1.12:0.02:1.2, 1.25:0.05:1.6, ...
    1.7:0.1:2, 2.2, 2.4, 2.8, ...
    3.6, 5, 7.5, 10];

% Delta Tao
dt = T/steps;
% Number of Asset Prices Grid = Q
Q = length(S);

% Now if we want to increase size of S grid
if NodeM ~= 1
    for i = 2:NodeM
        S_old = S;
        Q = length(S_old);
        S_new = zeros(1,2*Q-1);
        S_new(1:2:end) = S_old;
        S_new(2:2:end-1) = (S_old(1:end-1) + S_old(2:end))/2;
        S = S_new;
        Q = length(S);
    end
end

% Option Value at tao = T-T = 0, Default for put option
V = zeros(Q,1);
V(Q) = 0;
V(1:end-1) = max(S(1:end-1)-K1,0) - 2*max(S(1:end-1)-(K1+K2)/2,0)+...
    max(S(1:end-1)-K2,0);

% Now iterate through time steps to solve Vn
iteration = 0;
for i=1:steps
    [delta, gamma] = Greeks(S', V, method, 0, 0);
    M = uncertainM(gamma, S, r, sigmaMax, sigmaMin, cases); 
    M_new = M;
    % Now use Newton's Method to get Vn+1
    V_old = V;
    I = speye(Q);
    ratio = tol+1;
    while ratio > tol
        V_new = (I + (1-theta)*M_new*dt)\(I - theta*M*dt)*V;
        Top = abs(V_new - V_old);
        Bot = max(1, abs(V_new));
        ratio = max(Top./Bot);
        V_old = V_new;
        [delta, gamma] = Greeks(S', V_new, method, 0,0);
        M_new = uncertainM(gamma, S, r, sigmaMax, sigmaMin, cases);
        iteration = iteration + 1;
    end
    V = V_new;
end

% Now interpolate through V to find price based on K

Price = interp1(S,V,S0,'nearest');

if cases == 'w'
    [delta, gamma] = Greeks(S',V, method, 0);
end

end

function M = uncertainM(gamma, S, r, sigmaMax, sigmaMin, cases)

indGamma = (gamma <=0)'; % Index where Gamma < 0
% Calculate both Alpha and Beta based on sigmaMax and min
[alphaMax, betaMax] = alphaBeta(S, r, sigmaMax);
[alphaMin, betaMin] = alphaBeta(S, r, sigmaMin);
% Worst case long
if cases == 'w'
    Alpha = alphaMax.*indGamma + alphaMin.*(1-indGamma);
    Beta = betaMax.*indGamma + betaMin.*(1-indGamma);
else
    Alpha = alphaMin.*indGamma + alphaMax.*(1-indGamma);
    Beta = betaMin.*indGamma + betaMax.*(1-indGamma);

end

Q = length(S); 
% Now create the sparse matrix M with corrected Alpha and Beta
i = 2:Q-1;
M = sparse(i,i, (Alpha + Beta + r), Q, Q, 3*Q) + ...
    sparse(i,i-1, -Alpha, Q, Q, 3*Q) + ...
    sparse(i,i+1, -Beta, Q, Q, 3*Q);
M(1,1) = r;

end
